### =========================================================================
###
### IV. Occurrence cleaning with reference database ####
### 
### Please contact Kew (https://powo.science.kew.org/) for reference database
### Or you can use your own reference database
###
### The reference data should contain two parts:
### 1. Vector map (e.g. shapefiles)
### 2. data.frame with at least two columns indicating:
###   (1) species names
###   (2) names of distribution areas of this species, 
###       and the names (as well as its column name) should be the same with those in the related layers of the vector map.
###   e.g.
###     species      regions
###       sp1         area1
###       sp1         area2
###       sp1         area3
###       sp2         area1
###       sp2         area4
### 
### =========================================================================


### 1. Settings ####

### Please load reference datasets and map from Kew or your own database
refdata <- read.table("./Kew_data_world_plant/Plant_Data_Kew.txt", header=T)
refmap <- readOGR(dsn = "./Kew_data_world_plant/Kewmaps/level3.shp")


### Please load this function which corrected bugs in original CoordinateCleaner package.

source("./functions/clean_coordinates_debug.R")

### Set the working path and output folders for name correcting and coordinate cleaning 

wkpath <- ""

path.species <- "./2.3_cleaning_reference/"

### Load your occurrence data. 
### PLEASE ensure that there should include three columns: "name", "x", "y", 
### representing original species names, longitude, and latitude
spcleaned <- readRDS("./2.2_cleaning_cc/cleaning_cc_occurrences/Occurrence_cleaned.rds") # Example


#### 2. Preparation ####

### install packages globally
# set a local mirror server
options(repos=structure(c(CRAN="https://stat.ethz.ch/CRAN/")))

# Packages to load
packages <- c("rgdal", "sp")

# Install missing ones and load all packages
for (p in packages) {
  if(!library(package = p, logical.return = TRUE, character.only = TRUE)){
    install.packages(p)
    library(package = p, character.only = TRUE)
  } else {   
    library(package = p, character.only = TRUE) 
  }
}

### Set working path 
setwd(wkpath)
if (!dir.exists(path.species)){ #Create paths for cleaned species
  dir.create(path.species)}

### Prepare Kew dataset
refdata <- na.omit(refdata)
refdata_native <- refdata[refdata$Introduced == 0,]

refdata_family <- unique(refdata_native$Family)
refdata_native$name <- paste(refdata_native$Genus, refdata_native$Species)
refdata_name <- unique(refdata_native$name)


#### 3. Start cleaning ####

### Processing Kew Data ###

spfamily <- unique(occ.corrected$family)

spcleaned.df.all <- data.frame(stringsAsFactors = FALSE)
for (f in spfamily){
  
  print(paste("Working with:", spfamily[f]))
  
  #Create folders for cleaned species
  if (!dir.exists(file.path(path.species, spfamily[f]))){ 
    dir.create(file.path(path.species, spfamily[f]))}
  
  cleanrecord <- data.frame(stringsAsFactors = FALSE)
  
  if (spfamily[f] %in% refdata_family){
    
    #### =============== 1.1.1 Species in reference database =============== ####
    
    # Extract family from reference database
    refdata_family_table <- refdata_native[refdata_native$Family == spfamily[f],]
    
    # Extract family from cleaning database
    cleanfamily <- spcleaned[which(spcleaned$family == spfamily[f]), ]
    speciesname <- unique(cleanfamily$name)
    
    for (i in 1:length(speciesname)){
      
      print(paste("Working with", speciesname[i]))
      
      species <- cleanfamily[cleanfamily$name==speciesname[i], ]
      
      ### For species in the reference database
      if (speciesname[i] %in% refdata_name){
        
        # Check if the records fall in the Kew polygons
        coordinates(species) <- ~ x + y
        proj4string(species) <- proj4string(refmap)
        
        kewsprange <- refdata_family_table[refdata_family_table$name == speciesname[i], ]
        kewsprangemap <- refmap[which(refmap@data$LEVEL3_COD %in% unique(kewsprange$Area_code_L3)),]
        kewsprangemap2 <- gBuffer(kewsprangemap, byid=FALSE, width=2) # 2 degree
        
        sp_cleaned <- species@coords[names(na.omit(over(species, kewsprangemap2))),]
        
        # Cleaning summary
        if (class(sp_cleaned) =="numeric"){
          sp_cleaned <- as.data.frame(matrix(sp_cleaned,ncol=2))
          colnames(sp_cleaned) <- c("x","y")
          cleanrecord <- rbind(cleanrecord,
                               data.frame(species = as.character(speciesname[i]),
                                          original = nrow(species@coords),
                                          kept = nrow(sp_cleaned),
                                          deleted = nrow(species@coords)-nrow(sp_cleaned),
                                          delPercent = ((nrow(species@coords)-nrow(sp_cleaned))/nrow(species@coords)),
                                          status = "Checked",
                                          stringsAsFactors = FALSE),
                               stringsAsFactors = FALSE)
          
        } else { # for species with more than 1 record
          cleanrecord <- rbind(cleanrecord,
                               data.frame(species = as.character(speciesname[i]),
                                          original = nrow(species@coords),
                                          kept = nrow(sp_cleaned),
                                          deleted = nrow(species@coords)-nrow(sp_cleaned),
                                          delPercent = ((nrow(species@coords)-nrow(sp_cleaned))/nrow(species@coords)),
                                          status = "Checked",
                                          stringsAsFactors = FALSE),
                               stringsAsFactors = FALSE)
        }
        
        # Output cleaning results for each species
        write.csv(sp_cleaned, file.path(path.species, spfamily[f], paste0(spfamily[f], "_", speciesname[i],".csv")), row.names = FALSE)
        
        spcleaned.df <- as.data.frame(sp_cleaned)
        spcleaned.df$family <- spfamily[f]
        spcleaned.df$name <- speciesname[i]
        spcleaned.df.all <- rbind(spcleaned.df.all, spcleaned.df)
        
      }
      
    } else {
      
      ####=============== 1.1.2 Species not in reference database ===============####
      
      write.csv(species, file.path(path.species, spfamily[f],paste0("NotCleaned_", spfamily[f], "_", speciesname[i],".csv")), row.names = FALSE) #Names in my data but not in Kew
      
      spcleaned.df <- species[ ,c("x", "y", "family", "name")]
      spcleaned.df.all <- rbind(spcleaned.df.all, spcleaned.df)
      
      # Cleaning summary  
      if (class(species) =="numeric"){
        cleanrecord <- rbind(cleanrecord,
                             data.frame(species = as.character(speciesname[i]),
                                        original = 1,
                                        kept = 1,
                                        deleted = 0,
                                        delPercent = 0,
                                        status = "Not_Checked",
                                        stringsAsFactors = FALSE), 
                             stringsAsFactors = FALSE)
        
      } else {
        
        cleanrecord <- rbind(cleanrecord,
                             data.frame(species = as.character(speciesname[i]),
                                        original = nrow(species),
                                        kept = nrow(species),
                                        deleted = 0,
                                        delPercent = 0,
                                        status = "Not_Checked",
                                        stringsAsFactors = FALSE), 
                             stringsAsFactors = FALSE)
        
      }
      
      saveRDS(cleanrecord, file.path(path.species, "CleanRecord_byReference.rds"), row.names = FALSE)
      
    }
    
    print(paste("End of family:", spfamily[f]))
    
  } else { # Family not in reference database
    
    print(paste("No reference database for family:", spfamily[f]))
    
    ####=============== 1.2 Families and species not in Kew ===============####
    
    # Extract family from cleaning database
    cleanfamily <- spcleaned[which(spcleaned$family == spfamily[f]), ]
    speciesname <- unique(cleanfamily$name)
    
    cleanrecord <- data.frame(stringsAsFactors = FALSE)
    for (i in 1:length(speciesname)){ # Species from family not in reference database
      
      print(paste("Working with", speciesname[i]))
      
      species <- cleanfamily[cleanfamily$name==speciesname[i], ]
      write.csv(species, file.path(path.species, spfamily[f], paste0("NotCleaned_", spfamily[f], "_", speciesname[i],".csv")), row.names = FALSE) 
      
      spcleaned.df <- species[ ,c("x", "y", "family", "name")]
      spcleaned.df.all <- rbind(spcleaned.df.all, spcleaned.df)
      
      # Cleaning summary
      
      if (class(species) =="numeric"){
        cleanrecord <- rbind(cleanrecord,
                             data.frame(species = as.character(speciesname[i]),
                                        original = 1,
                                        kept = 1,
                                        deleted = 0,
                                        delPercent = 0,
                                        status = "Not_Checked",
                                        stringsAsFactors = FALSE), 
                             stringsAsFactors = FALSE)
        
      } else {
        
        cleanrecord <- rbind(cleanrecord,
                             data.frame(species = as.character(speciesname[i]),
                                        original = nrow(species),
                                        kept = nrow(species),
                                        deleted = 0,
                                        delPercent = 0,
                                        status = "Not_Checked",
                                        stringsAsFactors = FALSE), 
                             stringsAsFactors = FALSE)
      }
      
    } # Species from family not in reference database
    
  } # Family not in reference database
  
  saveRDS(cleanrecord, file.path(path.species, "CleanRecord_byReference.rds"), row.names = FALSE)
  
}

saveRDS(spcleaned.df.all, file.path(path.species, "Occurrence_cleaned_cc_kew.rds"))


tmp <- tempfile()
do.call(file.remove, list(list.files("tmp", full.names = TRUE)))
